Predicting Milling Stability Based on Composite Cotes-Based and Simpson’s 3/8-Based Methods

Avoiding chatter in milling processes is critical for obtaining machined parts with high surface quality. In this paper, we propose two methods for predicting the milling stability based on the composite Cotes and Simpson’s 3/8 formulas. First, a time-delay differential equation is established, wherein the regenerative effects are considered. Subsequently, it is discretized into a series of integral equations. Based on these integral equations, a transition matrix is determined using the composite Cotes formula. Finally, the system stability is analyzed according to the Floquet theory to obtain the milling stability lobe diagrams. The simulation results demonstrate that for the single degree of freedom (single-DOF) model, the convergence speed of the composite Cotes-based method is higher than that of the semi-discrete method and the Simpson’s equation method. In addition, the composite Cotes-based method demonstrates high computational efficiency. Moreover, to further improve the convergence speed, a second method based on the Simpson’s 3/8 formula is proposed. The simulation results show that the Simpson’s 3/8-based method has the fastest convergence speed when the radial immersion ratio is large; for the two degrees of freedom (two-DOF) model, it performs better in terms of calculation accuracy and efficiency.


Introduction
The three main types of vibration that occur during high-speed milling are free vibration, excited vibration, and self-excited vibration, and the regenerative chatter in self-excited vibration is the main cause of instability in the machining process [1]. Because this chatter typically leads to poor surface quality of the machined parts, aggravated tool wear, and even reduced service life of machine tools, it must be avoided to solve or overcome these problems [2]. In the analysis of chatter, dynamic milling processes that consider the delay effect are generally described using delay differential equations with respect to the timeperiodic coefficients [3,4]. The chatter stability based on these differential equations is an important index for achieving high-performance machining [5].
In recent decades, experimental, analytical, and numerical methods have been studied to obtain stability lobe diagrams, which can be used to determine exact values avoiding chatter. Wu et al. [6] utilized the Lyapunov index to measure whether chatter occurred; however, they could only predict the system stability under specific parameters. Davies et al. [7] used the time-domain calculation method to predict unstable regions in the stability lobe diagrams, which was only applicable to small radial depths of cut. Altintas and Budak et al. [8] first presented a zero-order approximation method to quickly obtain the stability lobe diagrams, where the real and imaginary parts of the feature are used to determine the cutting parameters. It is worth mentioning that this method was only applicable to high radial depths of cut. Bayly et al. [9] employed the time-finite element method to calculate tool motion and utilized the weighted residual method to determine the Floquet transition matrix. Based on the delay differential equations, Insperger et al. [10] proposed the semidiscretization method (SDM) and the first-order semi-discretization method (1 st SDM), which are widely used as the basis for evaluating other stability analysis methods in the time domain. The SDM and 1 st SDM only discretized delay states and periodic coefficients. Based on direct integration, Ding et al. [11] presented a full-discretization method (FDM) to predict milling stability. Unlike the SDM and 1 st SDM, the FDM carried out the linear interpolation on the state terms to improve computational efficiency. Further, Insperger [12] compared the FDM with the SDM and 1 st SDM in terms of the convergence speed and computational efficiency. It was found that the convergence speed of the 1 st SDM was faster than those of the SDM and FDM, while the computational efficiency of the FDM was higher than those of the SDM and 1 st SDM. After the comparisons, some FDM-based methods were used to predict the milling stability. In addition to the above methods, the numerical integration methods are widely utilized to obtain the stability lobe diagrams. Ding et al. [13] approximated the integral terms with the Newton-Cotes and Gauss integral formulas to predict the milling stability. Lu et al. [14] approximated the solution of the delay differential equation with the direct integration technique. Qin et al. [15] approximated the state terms of the delay differential equation with the Chebyshev wavelets of the second kind. Moreover, Liu et al. [16] combined the numerical solution of the delay differential equation with the Simpson formula to predict milling stability (SEM), which was only suitable for large radial immersion. However, the aforementioned numerical methods cannot simultaneously guarantee fast convergence speed and high computational efficiency.
To this end, we present two numerical analysis methods, namely the composite Cotesbased method (CCM) and Simpson's 3/8-based method (S38M), in this paper to improve the convergence speed and computational efficiency at the same time. The remainder of this paper is organized as follows. Section 2 presents the mathematical model of system motion with two degrees of freedom (DOF) and the state-space. The composite Cortesbased method and its simulation and analysis for the single-DOF model are presented in Section 3. Subsequently, the Simpson's 3/8-based method and its simulation and analysis for the single-DOF and two-DOF models are presented in Section 4. Finally, the conclusions are stated in Section 5.

Mathematical Model
Taking down-milling as an example, the dynamic system of end-milling with two degrees of freedom is shown in Figure 1, where the workpiece is assumed to be rigid, and the milling cutter is assumed to be evenly distributed. Note that the helix angle of the milling cutter is not taken into account. Considering the regenerative effect, the mathematical model of system motion can be expressed as a second-order differential equation [17] as follows: where q(t), . q(t), and .. q(t) represent the displacement vector, the first derivative of q(t) w.r.t t, and the second derivative of q(t) w.r.t t, respectively. M, C, K, and a p represent the mass matrix, damping matrix, stiffness matrix, and depth of cut, respectively. T = 60/(Nv) represents the delay period of the delay differential equations, where N represents the cutter tooth number, and v represents the spindle speed. The radial cutting force coefficient matrix K c (t) can be expressed as follows: where h xx (t) = N ∑ j=1 g φ j (t) sin φ j (t) K t cos φ j (t) + K n sin φ j (t) 3 of 23 h xy (t) = N ∑ j=1 g φ j (t) cos φ j (t) K t cos φ j (t) + K n sin φ j (t) (4) h yx (t) = N ∑ j=1 g φ j (t) sin φ j (t) −K t sin φ j (t) + K n cos φ j (t) (5) and h yy (t) = N ∑ j=1 g φ j (t) cos φ j (t) −K t sin φ j (t) + K n cos φ j (t) (6) where K t and K n represent the tangential and normal cutting force coefficients, respectively. φ j (t) represents the angular position of tooth j, as shown in Figure 1, and is defined as follows: where Kt and Kn represent the tangential and normal cutting force coefficients, respectively.
( ) j t φ represents the angular position of tooth j, as shown in Figure 1, and is defined as follows: is a piecewise function used to judge whether the tooth j is cutting or not and is defined as follows: where st φ and ex φ represent the start and exit cutting angles of the tooth j, respectively.
In the up-milling process,   g φ j (t) is a piecewise function used to judge whether the tooth j is cutting or not and is defined as follows:

Composite Cotes-based method
where φ st and φ ex represent the start and exit cutting angles of the tooth j, respectively. In the up-milling process, φ st = 0 and φ ex = arccos(1 − 2a/D); in the down-milling process, φ st = arccos(2a/D − 1) and φ ex = π, where a/D represents the radial immersion ratio.

Composite Cotes-Based Method
. Then, the state expression of (1) can be expressed as follows: .
where x(t) = Ax(t), then the following equation can be derived: (12) where t 0 represents the start time. It is worth noting that x(ξ) is unknown such that (12) is a Volterra integral equation of the second kind.

Numerical Algorithm
To solve (9) using the composite Cotes formula, the integral equation f (x) is assumed to be continuous in the interval [c, d]. Then, [c, d] is divided into m isometric intervals; that is, the span of each interval h is equal to (d − c)/m. Based on the isometric nodes u k = c + (k − 1)/h, where k = 1, 2, . . . , m + 1, we obtain the following: where C (m) k represents the Cotes coefficients [18], expressed as follows: It is worth noting that C (m) k only depends on m. When m = 4, (14) can be expressed as follows: where T consists of the durations of the free and excited vibrations [15]. The durations of the free vibration and excited vibration can be defined as t f and t c , respectively. For the free vibration, B(ξ) = 0, and (12) can be expressed as follows: Therefore, the state equation can be solved when t = t 0 + t f , and it can be expressed as follows: For the excited vibration, its time lies in [t 0 + t f , t 0 + T]. The interval is divided into n isometric intervals, and the span of each interval l is equal to (T − t f )/n. The discrete time can be expressed as t i = t 0 + t f + (i − 1)h, where i = 1, 2, . . . , n + 1. According to the numerical integral solution of the Volterra integral equation of the second kind [19], we obtain the following: If [t i , t i + 1 ] is divided into four isometric intervals, then t i , t i+1/4 , t i+1/2 , t i+3/4 , and t i+1 can be derived. Combining (13) and (15), the integral equation can be expressed as follows: Upon substituting (19) into (18), we can see that x(t i+1/4 ), x(t i+1/2 ), and x(t i+3/4 ) cannot be solved directly. The barycentric Lagrange interpolation method [20] is introduced here to approximate x(t i+1/4 ), x(t i+1/2 ), and x(t i+3/4 ) and obtain the following: Substituting (19)- (22) into (18), x(t i ) can be expressed as follows: It is worth noting that B(t n+2 ) cannot be calculated directly using (9) when i = n. A Newton Cotes formula [13] is introduced here to calculate x(t n+1 ) and is expressed as follows: Combining (17), (23), and (24), the transition matrix can be constructed as follows: where B i represents B(t i ), i = 1, 2, . . . , n + 1. According to (25)-(27), the dynamic mapping of the discrete system can be determined as follows: . . .
where I 1 is a (2n + 1) × (2n + 1) identity matrix. Therefore, the transition matrix of the dynamic system in a tool pass cycle is obtained as follows: It is worth noting that the chatter stability must be determined using the Floquet theory. If the modulus of any eigenvalue in Φ 1 exceeds 1, the system is unstable; if the moduli of all eigenvalues in Φ 1 are less than 1, the system is stable [21,22]. Therefore, the boundary curve between the unstable and stable regions in the lobe diagram can be used as the criterion for judging whether chatter occurs.

Simulation and Analysis
In this section, SDM [10] and SEM [16] are compared. For objective comparison, the CCM, SDM, and SEM share the same parameters and machining conditions. A benchmark example of the single-DOF milling model is utilized to validate and analyze the CCM. The single-DOF milling mathematical model is expressed as follows [23]: ..
where ζ denotes the relative damping, ω n denotes the angular natural frequency, and h(t) denotes the specific cutting force coefficient. The state-space for single-DOF milling mathematical model is expressed as follows: . where The parameters in (32) and (33) are defined as follows: f n = 922 Hz, ξ = 0.011, m t = 0.3993 kg, K t = 6 × 10 8 N/m 2 and K t = 2 × 10 8 N/m 2 . The adopted machining condition is down-milling. The radial immersion ratio a/D is set as 1 to avoid intermittent milling. The spindle speed is set as 5000 rpm (v = 5000 rpm). The depths of cut are set as 0.001 m, 0.0005 mm, and 0.0002 mm, respectively. All programs in this study are executed in MATLAB R2019a and run on a personal computer (AMD Ryzen 5 5600H; CPU 4.0 GHz, 16 GB). The maximum modulus of the eigenvalues of Φ is labelled as |λ|, and the maximum modulus of the eigenvalues of Φ with SDM is labelled as |λ 0 |. In this study, |λ 0 | is treated as the exact value.
The convergence rate comparisons of the CCM, SDM, and SEM are shown in Figure 2. It is seen from the figure that when n is small, the convergence rate of the CCM is significantly higher than those of the SDM and SEM regardless of the depth of cut. As n increases, the convergence rates of the three methods approach gradually, and the convergence rates of the CCM and SEM are higher than that of SDM. The figure (a) is enlarged to observe the differences better. It is worth noting that the convergence rate of the CCM is significantly higher than those of the SDM and SEM when the discrete number changes from 75 to 100. The results demonstrate that the convergence rate of the CCM outweighs those of the SDM and SEM. The stability lobe diagrams can be used to compare the calculation accuracy among the CCM, SDM, and SEM. As a reference, the discrete number is set as 500. The other parameters are set as follows: the discrete numbers are set as 25, 40, and 55; the radial immersion ratios are set as 0.05, 0.5, and 1; the spindle speed varies from 0.5×10 4 to 2.5×10 4 rpm, and 200 equally distributed sampling points are selected within this range; the depth of cut varies from 0 to 0.01 m, and 100 equally distributed sampling points are selected within this range.
In the single-DOF system, the stability lobe diagrams obtained with the CCM, SDM, and SEM when a/D = 1 are given in Table 1. It can be seen that the stability lobe diagrams of the CCM are closer to the reference ones than the SDM and SEM. Note that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25. Overall, the CCM is more accurate than the SDM and SEM. To further compare the calculation accuracy, two indicators are introduced, i.e., the arithmetic mean of relative error (AMRE) and mean squared error (MSE) [2]. The AMRE and MSE are defined as follows [2]:  The stability lobe diagrams can be used to compare the calculation accuracy among the CCM, SDM, and SEM. As a reference, the discrete number is set as 500. The other parameters are set as follows: the discrete numbers are set as 25, 40, and 55; the radial immersion ratios are set as 0.05, 0.5, and 1; the spindle speed varies from 0.5 × 10 4 to 2.5 × 10 4 rpm, and 200 equally distributed sampling points are selected within this range; the depth of cut varies from 0 to 0.01 m, and 100 equally distributed sampling points are selected within this range.
In the single-DOF system, the stability lobe diagrams obtained with the CCM, SDM, and SEM when a/D = 1 are given in Table 1. It can be seen that the stability lobe diagrams of the CCM are closer to the reference ones than the SDM and SEM. Note that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25. Overall, the CCM is more accurate than the SDM and SEM. To further compare the calculation accuracy, two indicators are introduced, i.e., the arithmetic mean of relative error (AMRE) and mean squared error (MSE) [2]. The AMRE and MSE are defined as follows [2]: where a i and a i0 denote the predicted axial depth and reference axial depth, respectively. n denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the AMRE and MSE are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the AMRE and MSE of the SEM cannot be attained. Take n = 40 as an example. The AMRE obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The MSE obtained with the CCM, SDM, and SEM is 2.62 × 10 −8 , 4.24 × 10 −8 , and 1.07 × 10 −7 , respectively. The results show that the CCM is closer to the reference values than the SDM and SEM.
where ai and ai0 denote the predicted axial depth and reference axial depth, respectively. n denotes the discrete number. The discrete numbers are set as 25, 30, 35, 40, 45, 50, and 55. The variations of the AMRE and MSE are depicted in Figure 3. Due to the fact that the severe distortions appear in the stability lobe diagrams obtained with the SEM when the discrete number is 25, the AMRE and MSE of the SEM cannot be attained. Take n = 40 as an example. The AMRE obtained with the CCM, SDM, and SEM is 0.041, 0.076, and 0.154, respectively. The MSE obtained with the CCM, SDM, and SEM is 2.62 × 10 −8 , 4.24 × 10 −8 , and 1.07× 10 −7 , respectively. The results show that the CCM is closer to the reference values than the SDM and SEM.                  tively. Furthermore, the stability lobe diagrams, AMRE and MSE, as well as the computational time when a/D = 0.05 and a/D = 0.5 are given in Tables 2 and 3 and Figures 5 and 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant advantages in the computational efficiency.     Tables 2 and 3 and Figures 5 and 6. It could be noted that the stability lobe diagrams obtained with the CCM are closer to the reference ones than those obtained with the SDM and SEM, and the CCM has the highest calculation accuracy. The results also demonstrate that the CCM has significant advantages in the computational efficiency.
tional time when a/D = 0.05 and a/D = 0.5 are given in Tables 2 and 3 an It could be noted that the stability lobe diagrams obtained with the CC reference ones than those obtained with the SDM and SEM, and the CC calculation accuracy. The results also demonstrate that the CCM h vantages in the computational efficiency.

State-Space Expression
To match with S38M, a new state-space expression is defined as follows [17]: .

x(t) = Ax(t) + a p B(t)[x(t) − x(t − T)]
where and Based on the numerical integration method of the Volterra integral equation of the second kind [19], the state-space equation can be deduced as follows: where t st represents the start time and defaults to t 1 . When t = t 3 , we obtain the following:

Numerical Algorithm
The Lagrange polynomial denoted by P L (x) [16,24] can be used to approximate f (x) as follows: where It is worth noting that x i , x i+1 , and x i+2 only exist when L = 2, and then, the Simpson's 3/8 formula (41) can be derived using P 2 (x).
x i+2 Based on (41), (38) can be rewritten as follows: (42) which can be expressed as follows: As B(t i+2 ) cannot be solved directly when i = n, the classical numerical integration method is introduced to determine x(t n+1 ) as follows: Combining (17), (43), and (44), the transition matrix can be constructed as follows: According to (45)-(48), the dynamic mapping of the discrete system can be determined as follows: Therefore, the transition matrix of the dynamic system in a tool pass cycle is obtained as follows:

Single-DOF Milling Model
A benchmark example of the single-DOF milling model is utilized to validate and analyze S38M. Its state-space expression is introduced as (31), where To analyze the convergence rate of S38M, the radial immersion ratio a/D is set as 1 to avoid intermittent milling. The depths of cut are set as 0.0003 m, 0.0008 m, and 0.0015 m. Figure 7 shows the convergence rate comparisons of the S38M, CCM, SDM, and SEM when the spindle speed is 8000 rpm. It is worth noting that the S38M exhibits the highest convergence rate compared with other chatter analysis methods. According to the enlarged figures on the right, the S38M has a good convergence rate even when the discrete number is large. Additionally, the results show that the convergence rate of the SEM fluctuates evidently, and the SEM is not as stable as the other methods. To further analyze the convergence rate of the S38M, the spindle speed is raised to 10,000 rpm, and the other parameters remain unchanged. Figure 8 shows the resultant convergence rates. It can be seen that these results are consistent with the previous results.
The low-immersion condition can be utilized to verify the stability of the convergence rate [25]. The radial immersion ratio is set as 0.05; that is, a/D = 0.05. The depths of cut are set as 0.0001 m, 0.0002 m, and 0.0003 m. The spindle speed was set to 5000 rpm. Figure 9 shows the convergence rate comparisons of S38M, CCM, SDM, and SEM under low immersion. The convergence rates arranged from high to low are in the following order: S38M, CCM, SDM, and SEM. Moreover, it is seen that the convergence rate of SEM fluctuated remarkably.  The low-immersion condition can be utilized to verify the stability of the convergence rate [25]. The radial immersion ratio is set as 0.05; that is, a/D = 0.05. The depths of cut are set as 0.0001 m, 0.0002 m, and 0.0003 m. The spindle speed was set to 5000 rpm. Figure 9 shows the convergence rate comparisons of S38M, CCM, SDM, and SEM under low immersion. The convergence rates arranged from high to low are in the following order: S38M, CCM, SDM, and SEM. Moreover, it is seen that the convergence rate of SEM fluctuated remarkably. The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.
As described above, the fluctuations appear in Figure 9. It is found that for the S38M, the convergence rate is often a local maximum when the discrete number is an odd one; the convergence rate is often a local minimum when the discrete number is an even one. To further study this problem, the discrete numbers are set as 30, 40, 50, 60, 70, and 80. The corresponding AMRE and MSE are plotted in Figure 12. The results show that the S38M attains a better computational accuracy when the discrete number is an even one.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.  The parameters in Section 3.2 are adopted for the S38M. Table 4 lists the computation time of the S38M. It is seen from the table that the S38M has a high calculation accuracy when a/D = 1. However, as a/D decreases, the calculation accuracy of the S38M is not as accurate as that of the CCM. The reason lies in that the interpolation of the S38M, which only uses three points, cannot reach a high calculation accuracy, while the interpolation of the CCM uses five points to interpolate, which made it calculate the smaller interval. The AMRE and MSE are shown in Figure 10. Notice that the calculation accuracy of the S38M and CCM is higher than those of the SDM and SEM when a/D = 0.05, and the calculation accuracy of the S38M is higher than those of the CCM, SDM, and SEM when a/D = 1. The reason for the former lies in that some large fluctuations exist in the S38M. The computational time of the above methods is compared in Figure 11. It is observed from the figure that compared with the SDM and SEM, the computational time can be significantly decreased with the S38M and CCM. Although the computational time of the S38M is almost equal to that of the CCM, the S38M has a higher convergence rate when the radial immersion ratio was large.    As described above, the fluctuations appear in Figure 9. It is found that for the S38M, the convergence rate is often a local maximum when the discrete number is an odd one; the convergence rate is often a local minimum when the discrete number is an even one. To further study this problem, the discrete numbers are set as 30, 40, 50, 60, 70, and 80. The corresponding AMRE and MSE are plotted in Figure 12. The results show that the S38M attains a better computational accuracy when the discrete number is an even one.

Two-DOF milling model
The two-DOF milling mathematical model is expressed as follows [17]:

t T a h t h t y t y t T
ζω ω ζω ω

Two-DOF Milling Model
The two-DOF milling mathematical model is expressed as follows [17]: x(t) ..
where the tool is symmetrical by default. m t , ζ, and ω n in the x-direction are identical to those in the y-direction. h xx (t), h xy (t), h yx (t), and h yy (t) are the same as those defined in where and The following deduction is the same as that in Section 4.1. The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

S38M SDM
a/D = 0.05 The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

Conclusions
The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

Conclusions
a/D = 0.5 The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

Conclusions
The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.  The following deduction is the same as that in Section 4.1.
The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

Conclusions
The following deduction is the same as that in Section 4.1.
The system parameters of the two-DOF milling model are identical to those of the single-DOF model. The radial immersion ratios are set as 0.05, 0.5, and 1. The figures shown in Table 5 are obtained over a 100×50-sized grid of the system parameters, that is, the spindle speed and the depth of cut. The thin, black, solid line denotes the reference stability lobe diagrams obtained using the SDM when the discrete number is 200 [26]. We observe that the accuracy of the stability lobe diagrams of the S38M is better than that of SDM. Additionally, the computation time is remarkably saved using the S38M.

Conclusions
In this study, we focused on the prediction of milling stability using composite Cotesbased and Simpson's 3/8-based methods. First, the composite Cotes-based method is proposed for preventing chatter in milling processes. The three steps for obtaining the milling stability lobe diagrams are as follows: (1) establish the time-delay differential equation while considering the regenerative effects; (2) determine the transition matrix using the integral equations; (3) analyze the system stability according to the Floquet theory. To measure the calculation accuracy, the AMRE and MSE are adopted in our work. The results demonstrate that for the single-DOF model, when the discrete number is 40, the AMRE and MSE can be respectively reduced from 0.076 to 0.041 and from 4.24 × 10 −8 to 2.62 × 10 −8 , and the calculation time can be reduced from 55.63s to 20.21s by using the proposed composite Cotes-based method. In addition, the Simpson's 3/8-based method is proposed to further improve the calculation efficiency. The results demonstrate that for the single-DOF model, the proposed Simpson's 3/8-based method significantly improves the convergence speed while sharing the same computation time with the composite Cotes-based method when the radial immersion ratio is large; for the two-DOF model, the accuracy of the stability lobe diagrams of the Simpson's 3/8-based method are better than those of the SDM, and the computation time is remarkably saved using the Simpson's 3/8-based method.